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The free surface lattice Boltzmann method (FSLBM) is a combination of the hydrodynamic 
lattice Boltzmann method (LBM) with a volume of fluid (VOF) interface capturing technique for 
the simulation of incompressible free surface flows. Capillary effects are modeled by extracting the 
curvature of the interface from the VOF indicator function and imposing a pressure jump at the free 
boundary. However, obtaining accurate curvature estimates from a VOF description can introduce 
signihcant errors. This article reports numerical results for three different surface tension models 
in standard test cases, and compares the according errors in the velocity field (spurious currents). 
Furthermore, the FSLBM is shown to be suited to simulate wetting effects at solid boundaries. To 
this end, a new method is developed to represent wetting boundary conditions in a least squares 
curvature reconstruction technique. The main limitations of the current FSLBM are analyzed and 
are found to be caused by its simplified advection scheme. Possible improvements are suggested. 
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I. INTRODUCTION 

The free surface lattice Boltzmann method (FSLBM) 
[T] is a numerical method for the simulation of free sur¬ 
face flows combining a volume of fluid (VOF) approach 
[HU for interface tracking with the lattice Boltzmann 
method (LBM) [SHS] for hydrodynamics. We use the 
same definition of free surface flow as described in 
where it denotes a single phase flow problem contain¬ 
ing free boundaries instead of a two phase flow problem. 
VOF methods follow the notion of a sharp interface repre¬ 
sentation, i.e., assuming hydrodynamic equations for the 
bulk of the flow and modeling the interface by boundary 
conditions or by a jump of flow parameters in one-fluid 
approaches for two-phase flows. This is in contrast to cur¬ 
rently popular lattice Boltzmann multiphase approaches 
[61IIIIII2] (e.g., color gradient model [13], Shan-Chen 
model [13], free energy model m). that are based on a 
diffusive interface assumption. In the FSLBM, the LBM 
is used only to approximate the incompressible Navier- 
Stokes equations for the liquid phase. With sharp in¬ 
terface simulation techniques capillary effects need to be 
modeled explicitly in addition to the interface tracking. 
Altogether, there are three components that the total 
accuracy of the method depends on: the hydrodynamic 
solver (LBM), the interface tracking (advection of indi¬ 
cator function), and the surface tension model. While 
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the accuracy of the LBM is well-understood [TMI^ , only 
very few results have been reported addressing the accu¬ 
racy of the FSLBM’s advection scheme or its validity to 
simulate surface tension. Nevertheless, the approach has 
found numerous applications in surface tension driven 
complex flow scenarios [TM^ and even high Reynolds 
number flows |24l I25j. Hence, in the present paper we 
evaluate the accuracy of the FSLBM in surface tension 
driven flows, and also discuss existing limitations due to 
the original advection scheme. While conventional VOF 
implementations rely on geometric reconstruction to ap¬ 
proximate the advection equation of the indicator func¬ 
tion 0111: the FSLBM instead exploits the specific na¬ 
ture of the LBM. However, the present results in this 
work indicate that though this simplification reduces the 
algorithmic complexity significantly, it also comes with a 
comparably low accuracy. 

The simulation of surface tension involves the extrac¬ 
tion of curvature information from the interface defined 
through the fill level function [51 E51B5] . The VOF indica¬ 
tor function is non-smooth by definition, making the esti¬ 
mation of the interface curvature a non-trivial task. Nu¬ 
merical errors in the determination of interface stresses 
lead to the effect of undesired spurious currents that can 
invalidate or destabilize the computed solutions. We re¬ 
mark, that, if the surface stress is included as force term 
in the momentum equation as in |26j , then the incompat¬ 
ible discretization of pressure and indicator function gra¬ 
dients becomes an additional source for spurious currents 
[29] . The present model does not have this problem, be¬ 
cause surface tension is imposed as a pressure jump at the 
boundary instead. Anyway, errors in the curvature esti- 
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mation are reported as problematic - primarily because 
the grid convergence is poor unless more sophisticated 
approaches are employed. Modern VOF codes therefore 
rely on higher order curvature estimation techniques to 
suppress this error inisn]. However, these techniques re¬ 
quire the use of larger stencils. This may be undesirable 
or is even impossible if only next-neighbor information 
is available, as e.g. in some parallel computing environ¬ 
ments. Hence, in this paper we review and present three 
different techniques to approximate the interface capil¬ 
lary tension from the fill level function data of a free 
surface code that are using local 3x3x3 neighborhoods 
only. 

The first method is an adaption of the classic finite 
difference (FD) model by Brackbill et al. [55] that uses 
finite difference computations to obtain the interfacial 
curvature. As reported previously for this method, the 
magnitude of the spurious currents makes it impossible 
to simulate small capillary number flows correctly. No¬ 
tice that in |26) . the surface tension appears as a force 
term in the momentum equation, and further that this 
force is smoothed out over several grid cells. Hence this 
model is also called continuous surface force model in the 
literature. In the FSLBM presented here, surface tension 
is included in the boundary condition instead. 

The second model locally reconstructs the interface 
described by the VOF function as a continuous surface 
made up of triangles. From this triangular reconstruction 
(TR) the curvature information can be extracted. This 
method has first been presented in |31j . and turns out 
to effectively reduce the error due to spurious currents, 
because the predicted curvature values are much more 
accurate. 

The third method also involves a local reconstruction 
of the interface geometry. Based on a least squares ap¬ 
proach a parabolic approximation to the interface is con¬ 
structed in a local neighborhood of cells. From this least 
squares reconstruction (LSQR) one obtains curvature es¬ 
timates with a high accuracy. A similar approach has also 
been described in [35] as a fallback solution for a higher 
order height-function technique. For the current paper, 
we have extended the approach to problems including ad¬ 
hesive boundary conditions. It is the only model in this 
study achieving a second order rate of convergence in the 
classic spherical bubble benchmark. 

We also evaluate the boundary conditions needed to 
simulate the wetting behavior of solid surfaces for the 
FD and LSQR model. The original work by Brackbill 
et al. [55] discusses adhesive boundary conditions, which 
we adopted to the FSLBM context for this paper. Since 
curvature estimation based on finite differences often in¬ 
troduces larger errors, smoothing and filtering techniques 
are typically used and have been extended to adhesive 
boundaries in [33] for the simulation of low Capillary 
number flows. Motivated also by the limited accuracy of 
the finite difference approximation of curvature by finite 
differences, m and [31] switched to the height function 
technique [3^ |35] for the simulation of surface wetting. 


Some other approaches use level sets instead of the VOF 
method to represent the interface [37], or directly use 
non-Eulerian techniques (e.g. [35]). In this paper, we 
present the results obtained from the FSLBM in several 
simulations of wetting surfaces while studying alternative 
surface tension models. This allows a comparison be¬ 
tween the FD approach to surface tension and the newly 
developed LSQR method. 

II. NUMERICAL METHOD 

Throughout this paper we assume a three-dimensional 
D3Q19 - lattice Boltzmann model [S] [39] with Q = 19 
lattice velocity vectors Cg with q = 0, — We denote 

the LBM data (discrete particle distribution function - 
PDF) by f = (/o, /i, /q-i)- The LBM data is defined 

for every node within the liquid subdomain H(f), and 
follows the evolution equation 

fq{x + Cgff + 1) = fg{x,t), (la) 

f'(x,t) = f(x,t)-t-C(f(x,t)). (lb) 

where C is a collision operator, and f' is the post-collision 
distribution function. For the present paper, we have 
used the hydrodynamic two relaxation time (TRT) col¬ 
lision operator [ID]. The TRT operator has two eigen¬ 
values Ae, Ao, controlling the relaxation of the even and 
odd parts of the PDF, respectively. Similar to the pop¬ 
ular LBGK model [39] the relaxation rate r = — 1/Ae 
controls the kinematic viscosity ly = Cg(T — 1/2), with 
Cs = I/'s/s for the present D3Q19 model. The second 
parameter Ao is chosen to minimize the error at straight 
axis-aligned walls. The macroscopic variables of pressure 
P(x, t) = c^p(x, t) and velocity u(x, t) are defined via 
moments of the distribution function, 

Q-i 

Q-i 

pu=Y^ Cgfg, (2b) 

9=0 

at the respective node x. 

The FSLBM has first been described in [T]. To track 
the interface position in simulations, a VOF approach 
introduces a fill level function ip following an advection 
equation over time. The function ip is defined for each 
finite volume (or cell) as the volume fraction filled with 
liquid. Only the flow within the liquid subdomain H(t), 
consisting of all cells x with positive fill level (/?(x) > 0, is 
simulated by means of the LBM. The remaining cells, 
referred to as gas cells, become temporarily inactive. 
Within the liquid subdomain H(t), we distinguish be¬ 
tween liquid and interface cells. To count as interface, 
the cell must have both liquid and gas cells in the direct 
neighborhood defined by the lattice model. Only the in¬ 
terface cells are allowed to have a fill level between zero 
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and one. The interface cells are also used to compute 
the advection of the free surface in terms of the fill level 
function, and to impose a free surface boundary condi¬ 
tion on the LBM. In comparison to other VOF methods, 
where the advection of the fill level function is a non¬ 
trivial problem that needs to be solved in addition to the 
hydrodynamics, the FSLBM exploits the nature of the 
lattice Boltzmann equation, Eq. (la I and Eq. , to 
update the fill levels directly [T], and sets 


+ 1 ) = 


q) ■ (/^(x -b c„ t) - /' (x, t)) 


p(x,t + 1) 


(3a) 


where the coefficient k{x,q) is 


{ 0 if X + Cg is gas, 

I [(^(x + Cq) + i^(x)] if X -b Cg is interface, 

1 if X + Cg is liquid, 

(3b) 

which means that there is no mass flux to gas cells. 

The free surface boundary condition for the LBM at 
the interface cells is 


/,j(x,t + 1) = fq'^ipb,Ub) + /|’(pb,Ub) - fq(TC,t), (4) 

for all directions q pointing towards gas cells x -b ^ 
Hereby, /' denotes the post-collision distribution 
oriented towards the gas phase and fq'^{pb, Ub) is the equi¬ 
librium distribution function [1]. It can be shown that 
Eq. @ approximates a free boundary condition with first 
order spatial accuracy m- The two boundary condi¬ 
tion parameters pb and Ub are the macroscopic pressure 
and velocity at the interface, respectively. The boundary 
value Pb is defined as 


Pb = ^iPg+Ps), (5) 

where Pg{t) is the static pressure at the free surface and 
Ps (x, t) is the Laplace pressure. The latter depends on 
the local curvature of the interface by 

Ps(x) = 2 (tk(x), (6) 


with a constant surface tension a. We remark, that in the 
original work of [1], Eq. Q is applied for all q, oriented 
outwards with respect to the interface, i.e., • n < 0, 

where n = V(/3 is a local normal vector to the interface 
pointing towards the liquid (cf. Sec. |II A ). However, we 
find that this approach leads to anisotropic artifacts in 
the free surface dynamics. Eigure documents one ac¬ 
cording example. Hence, for the present work, we take 
into account the interface orientation only at the contact 
line (corner nodes) where the boundary becomes inho¬ 
mogeneous. As shown in Fig. a contact line cell has 
both links to solid wall (off-boundary, S) nodes and to 


FIG. 1: Simulation of a droplet splashing onto a liquid 
film without surface tension. Reynolds number « 250. 



(a) Slice through domain (x = y). Left: Imposing the free 
boundary condition on all links c, • n < 0 does not 
reproduce the rim correctly. Right: Rim clearly visible with 
the present implementation. 



(b) Iso-contour of the tracked interface. Left: Imposing the 
free boundary condition on all links • n < 0 leads to 
incorrect crown formation. Right: Same simulation based on 
present implementation. 



Node types: 

L: liquid 
I: interface 
G: gas (off-boundary) 
S: solid (off-boundary) 


FIG. 2: Inhomogeneous boundary at an interface corner 
node with boundary-intersecting lattice directions 
indicated by arrows. The free boundary rule is imposed 
on the thick arrow direction due to its orientation with 
respect to the interface with normal n. 


gas (off-boundary, G) nodes. Links to solid wall nodes 
(off-boundary nodes S) are usually subject to a no-slip 
condition |42j . The present implementation replaces the 
no-slip condition with the free surface condition in cor¬ 
ner nodes, if the outgoing lattice direction has negative 
orientation with respect to the inward-oriented interface 
normal n projected into the solid wall. 

The remaining part of this section lays out the three 
different approaches considered in this study, to extract 
the local interface curvature k from the volume fraction 
ip. Notice that the notion of a lattice cell reflects the rep¬ 
resented cubic volume with the length of one grid spac¬ 
ing, and which is used to define the fill levels. However, 
the LBM data is more precisely located at lattice nodes 
which coincide with the centers of the cells. 
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A. Finite Difference Approximation (FD) 

As first proposed in [26], one way to approximate the 
curvature k of the boundary surface defined through the 
fill levels, is to compute a finite difference approximation 
to 

« = -(V-n), (7) 


neighborhood around each interface cell. The interface 
points are determined using a piecewise linear interface 
construction (PLIC) approach |5|. For an interface cell 
centered around x^, let P be the half space P — {x|(x — 
p) • n(xi) > 0}, where p = x^ + an(xi) is a point within 
the corresponding unit volume P(xi). Now, the interface 
point p is defined such that the cut-off volume VHP 
satisfies 


where n is the normalized gradient of the indicator func¬ 
tion. To obtain this gradient, we use central finite differ¬ 
ences to approximate 

n = V(/j, (8) 

which can be interpreted also as a local normal vector 
to the free surface. We use the finite difference approx¬ 
imation suggested in |4d] . However, the curvature com¬ 
putation according to Eq. Q means that second order 
derivatives of the non-smooth indicator function ip are 
approximated by finite differences, which inevitably in¬ 
troduces larger errors. Hence, much work has been pub¬ 
lished (cf. for instance [44] ) on effective ways to mollify 
the fill level information and smooth out surface tension 
in the “continuum surface tension” approach. We use 
the Ks - kernel with support radius e = 2.0 (cf. |44]b 
i.e., only the next neighbor information is included in 
the convolution. 

Wetting properties are included by directly specifying 
an ideal equilibrium contact angle Ogq for solid bound¬ 
aries. For obstacle cells with surface normal n^, the 
boundary condition at the solid wall is 

fi = fi^cos^e,-|-ntsin0eq, (9) 

where is a tangent vector to the wall and normal to 
the contact line |2S]. Since the wall position rarely co¬ 
incides with the lattice nodes, the boundary value ac¬ 
cording to Eq. (|^ is extrapolated to the obstacle node. 
Hereby, the vector ht is computed by projection of the 
interface normal at the fluid boundary cell x;, onto the 
wall. The boundary condition for the fill level in the ob¬ 
stacle cells (affecting the interface normal n(xh) in the 
boundary cells) is a reflection condition for the (^-values. 
Here, we generally compute the boundary value for the 
obstacle cells Xo ^ based on the neighboring inner 
nodes Xq -|- c, S using the formula 

(^(Xo) = ^ \h^ ■ Cq\p{Xo + Cq) / ^ |n^-Cq|, 

|n^-Cq|>a 

Xo+Cqef2 

( 10 ) 

with an apperture a = to achieve a smoothed 

reflection for the boundary values. 


B. Triangular Reconstruction (TR) based on 
piecewise linear interface construction (PLIC) 

In |2I] , a curvature computation is suggested based on 
a local triangulation of interface points in a 3 x 3 x 3 


vol{V n P) = p{xi). (11) 

We determine a iteratively, similar to [45] . with an er¬ 
ror bound of « 4 X 10“^^ (40 iterations in a bisection 
algorithm) assuming an exact surface normal. The inter¬ 
face point can be computed in one step together with the 
estimation of the surface normals. For the latter we em¬ 
ploy again a finite difference scheme to Eq. (|^, however, 
without any convolution step. Alternative, higher order 
PLIC algorithms are discussed in [351ISH] • 

Once the interface points are determined by the PLIC 
scheme, the algorithm described in [3T| is used to con¬ 
struct a local “triangle fan” from the interface points 
within the local neighborhood. Then, a variant of the al¬ 
gorithm described in |46j determines the curvature of this 
polygonal surface. Notice, that the described TR scheme 
as well as the curvature estimation by LSQR of Sec. |H Cj 
are based solely on the surface points and use the gradi¬ 
ent information represented by the interface normals only 
as far as it is needed to construct these interface points. 
The TR method can be extended to support adhesive 
boundary conditions. In m, a way to extend the local 
triangulations at solid boundaries to achieve an “artificial 
curvature” matching with a desired equilibrium contact 
angle, is described. However, the implementation of the 
geometric construction is difficult in three dimensions. 
Also, we found the results obtained from that method 
often not convincing, which motivated the least squares 
- based approach of the following section. 


C. Least Squares Reconstruction (LSQR) based on 
piecewise linear interface construction (PLIC) 


The third approach to include surface tension consists 
in reconstructing the interface as a quadratic function 
in each 3x3x3 neighborhood around an interface cell. 
It has previously been described in [33], there however 
without the inclusion of wall adhesion effects. Like the 
TR approach, it is base d on the PLIC of interface points 


described above in Sec. 


HB 


Let t„, t„ and n be a local 
tangential to the interface, 


orthonormal basis, i.e., 
and p the local interface point. Now, assume that i is 
indexing all the remaining interface cells ina3x3x3- 
neighborhood. The interface cell data (ni,Pi) is used to 
fit the model function 


/(u, v) = Av? -\- Bv^ -\- Cuv Hu Iv J (12a) 
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FIG. 3: Determination of contact point ^ for a 
contact line cell. The PLIC segment defined by 
interface point pi and interface normal is extended 
and intersected with the obstacle wall. 


with parameters (A, C, H, I, J), by minimizing the er¬ 
ror 


E = (12b) 

i 


Here, Ui = (pi - p) • t„, Vi = (pi - p) • t„, and /* = 
(Pi ~ P) ■ n- This yields a linear least squares problem 
that has to be solved locally. We obtain the best results, 
when fixing the constant parameter J = 0, i.e., accepting 
only solutions that interpolate the surface point p of the 
respective interface cell. We use the implementation from 
the LAPACK library [?5] based on QR decomposition of 
the corresponding system matrix. Once / is determined 
the curvature can be evaluated analytically, using 


A(1 P) + B{1 + H^) - 2CHI 
~ ' 


(13) 


The approach can be seen as a modified version of 
the PROST - scheme from [?S]. Both schemes fit a 
parabolic function in a local neighborhood around each 
interface cell. However, as a major difference, PROST 
fits f{x,y,z) directly to the fill levels minimizing the er- 
Si (^i(/) ~ of tfio cut-off volumes Vi that the 
iso-surface f{x, y, z) = 0.5 cuts out of the interface cell i. 
This makes the least squares problem non-linear and / 
has to be determined by iteratively computing the error 
and updating of coefficients. The scheme described here 
is computationally less expensive. 

To include the effect of boundary adhesion, we extend 
the method in the following way: If the local 3x3x3 
-neighborhood contains an obstacle cell, then for each 
contact-line cell (i.e., an interface cell that has an obsta¬ 
cle cell as neighbor) from the same neighborhood, one 
contact point is approximated with a contact normal rh 
defined according to Eq. In the contact point, we 
require 


V/(mc, Vc) = -{rriu, my)jmn, (14) 

where mu = my = m-t„, and m„ = m-n„, and Uc, 

Vy are the coordinates of the contact point in the locally 
defined tangential plane. If the neighboring interface cell 


FIG. 4: The bottom-left interface (I) cell in Fig. 4a has 
several gas (G) neighbors but no liquid neighbor(L). 
Hence, the configuration of Fig. is not supported by 
the present implementation. Shown in Fig. |4b| is a valid 
configuration since all interface cells have both liquid 
and gas neighbors. The minimum supported film 
thickness is therefore at least one (liquid) lattice cell. 



is the center of the current 3x3x3- neighborhood, we 
use Eq. as a constraint to the respe ctive optimization 

12b). Otherwise, the 


problem given by the Eqs. ( 

12a 

12b) 

condition is simply included 

in t 

le 0] 


:he optimization of the 
error, Eq. dT^ . To obtain the contact point, we con¬ 
struct the closest intersection of the interface segment of 
the contact line cell with the solid surface as in Fig. 


D. Current limitations 


The interface tracking of the present implementation 
defines interface cells as active lattice Boltzmann cells 
that have a D3Q19 neighborhood containing both gas 
and liquid cells. This definition turned out to impose a 
limitation when simulating thin liquid films on wetting 
surfaces and with contact angles below 45°. As shown in 
Fig.i the thickness of a liquid film on solid substrate has 
to be resolved at least by one liquid cell in height, for the 
method to work correctly. A film thickness smaller than 
I lattice unit is not supported. For strongly wetting sur¬ 
faces {9eg < 45°), we often observe anisotropic errors be¬ 
cause it is then problematic to impose the correct contact 
angles according t o the LSQR method (Sec. H C) and the 
TR method (Sec. HB): Depending on the approximated 
interface position represented by the fill level informa¬ 
tion, the computed curvature values then tend to oscillate 
and overshoot. In Fig. |4b[ for instance, the approxima¬ 
tion of the interface as a smooth surface through the two 
leftmost interface cells can be expected erroneous under 
the condition of an acute intersection angle {9eg < 45°) 
with the solid surface. It is important to notice that 
these restrictions are specific to the presented FSLBM 
algorithm, while the presented curvature reconstruction 
schemes can be applied in any VOF-context. 
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III. NUMERICAL RESULTS 


The numerical results presented in the following have 
been obtained with the waLBerla lattice Boltzmann 
framework that includes the FSLBM implementa¬ 
tion described in [3T]. 


A. Equilibrium Spherical Bubble 


The standard benchmark for surface tension models 
is a static equilibrium bubble. If the curvature estima¬ 
tion would return the exact value k = 1/R everywhere 
on the interface, the solution to the problem would be 
a perfectly vanishing velocity field. Due to the existing 
errors, however, spurious currents occur around the inter¬ 
face from regions with overestimated Laplace pressure to 
positions underestimating the value (cf. Fig. [5]. For 
non-sophisticated methods, the magnitude of the spuri¬ 
ous velocities can be related to the ratio cr//r, of surface 
tension and dynamic viscosity /i = pv, with a constant 
prefactor < 1. This means that the error is independent 
of the spatial resolution or converging very slowly, and 
thus poses a limitation in terms of applicability to prob¬ 
lems involving dominant capillary forces. See [IHHHISI] 
for a discussion of the problem in connection with the 
VOF method. 


Here we adopt the problem stated in [49] using non- 
dimensional values with respect to a cubic domain of unit 
length, containing a spherical bubble of radius R = 0.125 
centered at (0.5,0.5,0.5). We apply no-slip boundary 
conditions at top {z = 1) and bottom z = 0 sides of the 
domain and periodicity along all other directions. The 
viscosity of the liquid of density p = 4 is set to /r = 1, 
and the surface tension parameter is ct = 0.357. The di¬ 
mensionless Ohnesorge number {Oh = p/\/apR) corre¬ 
sponding to the given problem is Oh ~ 2.37. Notice that 
in |1H] the bubble is actually a second fluid of the same 
density and viscosity as the liquid (two-phase), while in 
our case the flow inside the gas bubble is not simulated 
but represented only in terms of a gas pressure value that 
is dynamically adjusted according to the changes of the 
gas volume upon interface advection (free surface flow 
with bubble model [5H|53]). We remark that the present 
test case is appropriate to evaluate the error in the cur¬ 
vature computation only. Since there is no flow velocity 
in this test case aside from the spurious currents, one ex¬ 
pects no significant error contribution from the advection 
of the indicator function or the LBM. We have therefore 
evaluated the test case first in the standard case with a 
static frame of reference (Sec. Ill A 1 ), and second in a 
moving frame of reference (Sec. [Ill A 2 ) with a constant 
uniform background velocity added to the flow. 



FIG. 5: Slice through a domain containing a single gas 
bubble. The color indicates the magnitude (lattice 
units) of the spurious currents due to errors in the 
curvature estimation with the FD-approach. 



FIG. 6: Temporal evolution of maximal and average 
(spurious) velocity for the different surface tension 
models at a fixed grid spacing 5^ = 1/96. After the 
decay of an initial shock, the magnitude of the spurious 
currents in the system can be evaluated. The FD model 
shows the largest errors and often leads to oscillative 
behaviour. Similar behavior is obtained for 6^ = 1/48, 
1/144, and 1/192. 
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FIG. 7: Dependency of spurious velocity (maximum and average) on grid spacing Sx for the three different models 
evaluated after 500, 000 time steps. The dotted and dashed lines represent first and second order slopes, respectively. 
For the FD scheme, the plot indicates a first order convergence for both maximum and average spurious flow speed. 
The errors of the FD scheme are orders of magnitude above those of the methods based on surface reconstruction. 


1. Static frame of reference 


We perform a resolution study varying the grid spac¬ 
ing between the values Sx = 1/48,1/96,1/144,1/192 
at a fixed time step of St = 10“"*, which yields the 
resolution-dependent lattice relaxation times r = 0.6728, 
1.191, 2.055 and 3.264. The surface tension parameter 
in lattice units varies accordingly and takes the values 
ctl/ 10-3 = 0.0987, 0.7896, 2.665 and 6.317. For ini¬ 
tialization of the fill levels at t = 0 with the spherical 
geometry, we employ a spatial subdivision technique, re¬ 
fining each discrete cell volume by a factor of 100 along 
each coordinate. The simulation then exhibits a series of 
pressure disturbances until the numerical equilibrium is 
reached. Fig. shows the development of the maximal 
and average flow velocity within the domain for 500,000 
time steps. After a certain number of time steps the 
shock wave is sufficiently decayed for both the TR and 
the LSQR method, such that both maximal and aver¬ 
age flow velocity within the domain become smaller than 
< 10“^° in magnitude. The FD approach, however does 
not converge and enters into an oscillating behavior in¬ 
stead. Here, the spurious currents are large enough to 
trigger changes in the layer of interface nodes. This also 
introduces sudden changes in the curvature computation, 
thus explaining the oscillations. Fig. shows the maxi¬ 
mal and average velocity within the domain for different 
spatial resolutions and various methods. The strength of 
the spurious currents obtained with the FD method are 
in accordance with the values reported in |49j . where a 
similar approach (CSF) is used for referencing. The cur¬ 
vature information obtained by geometric reconstruction 
(TR and LSQR methods) is much more accurate than 
the FD approximation, and reduces the spurious veloci¬ 


ties almost down to the order of machine precision. 

We also evaluate the accuracy of the curvature val¬ 
ues obtained for the numerical equilibrium, i.e., the state 
reached after 500,000 time steps. Fig. compares the 
error in curvature at various grid spacings for the three 
different methods. Only the plot for the LSQR method 
indicates a second order rate of convergence. However, 
it is clear that the occurrence of spurious currents is not 
due to constant over- or underestimation of curvature, 
but rather because of its variance with the node position. 
This becomes obvious, when comparing Fig. |8b| showing 
the standard deviation over all interface nodes in the final 
state to the resulting spurious currents of Fig. 

2. Moving frame of reference 

When dealing with dynamical problems numerical er¬ 
rors in the advection of the indicator function ip are of¬ 
ten critical. This holds in particular for the present test 
case, since any errors in the fill levels will introduce an 
additional error into the computed curvature values. To 
study effects of advection, we add a constant background 
velocity of Ug = (1,0,0)^ to the test case. This means 
that there is now a constant advection involved and the 
spurious currents appear as a deviation from the back¬ 
ground velocity. For Sx = 1/96, St = 2.5 x 10“® (lattice 
relaxation time r = 0.6728, Ux = Q.QQ2SSx/St), we run 
the simulation for T = L/uq^x time steps, i.e., the bub¬ 
ble traverses the periodic domain of length L exactly one 
time. This can be interpreted as a moving frame of ref¬ 
erence while, physically, the setup is equivalent to the 
static version of Sec. |HIA1| Measuring the curvature er¬ 
ror over time. Fig. [^exhibits a dramatic increase in error 
as compared to the static test case. The errors of the 




















FIG. 8: Comparison of curvature estimation for a stationary bubble after 500,000 time steps. For the FD-model, 
the included graphs are somewhat arbitrary, because the values oscillate over time, analog to the spurious currents 

in (cf. Fig.lol. 




grid spacing 

(a) The norm of the curvature error. Only the LSQR 
curvature error does converge in the test, with a rate of 
convergence in 0{S^). 


grid spacing Sx 

(b) Standard deviation of average curvature values. The FD 
scheme has the largest standard deviation. Consecutively, 
the FD scheme generates the largest spurious currents. 


reconstructive methods, LSQR and TR, are now signifi¬ 
cantly larger than the error of the FD model. A possible 
reason is that the FD model is more diffusive than the re¬ 
construction methods, and thus less sensitive to errors in 
the indicator function field. A grid study with space steps 
Sx = 1/48,1/96,1/144,1/192 and corresponding time 
steps (5t = 1 X 10-4,2.5 x 10-^ 1.11 x 10-^6.25 x 10"® 
(diffusive scaling) revealed that these errors do not con¬ 
verge with the grid spacing. Figure [T^ shows that the 
shape of the bubble after advection deviates notably from 
a true sphere. In accordance with the increased errors 
in curvature, the reconstruction based schemes show the 
most deviation. Also, in this case the spurious currents 
do no longer converge for either method. This indicates 
that there is an additional error that stems from the ad¬ 
vection of the indicator function ip. 


Since we are not aware of any numerical evaluation 
of the FSLBM advection scheme, we supplement the 
Laplace test with a convergence check of the indicator 
function values. To this end, a gas bubble of diam¬ 
eter d = lOda; is placed inside of a periodic compu¬ 
tational domain F with a prescribed uniform velocity 
Uq = (0.05,0, 0)’^(53;/i 5(. The LBM data is thus constant 
with f(x, t) = f®'^(pg/cg, Uo), where Pg is the constant 
reference pressure. This excludes any error contribution 
by the LBM or the surface tension modeling to the in¬ 
dicator function ip. From the prescribed LBM data we 
compute the advection of the bubble in terms of the in¬ 
dicator function according to Eqs. (3a I and (3b). After a 
time T = d/uo^x the bubble has moved a distance equal 
to its diameter. To evaluate the error, we use the and 
error norms, by comparing the fill levels at time T to 



FIG. 9: - error in curvature for the three different 

surface tension models for a moving frame of reference. 
Simulation of a spherical bubble during traversal a 
periodic domain with uniform velocity. At < = 0 the 
bubble is initialized to a nearly ideal (spherical) shape. 
Due to errors in the advection operator and the surface 
tension models the error increases and after t = 0.4 
oscillates around a fixed value. The FD scheme is less 
sensitive to errors in the fill levels, presumably because 
it is more diffusive and based on the mollified indicator 
function. 


STDEV{k) 
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FIG. 10: Comparison of the shapes of the bubbles after 
advection along a;-axis with surface tension. From left 
to right: FD, LSQR, TR. Visualization shows the layer 
of interface cells in the x-y-plane through the center of 
the bubbles at time 0.94r (150000 time steps). The 
initial radius of the bubble was 24Sx- 


O 


XI 


the initial configuration t = 0, i.e., 

Sxgr |y’(x + Tuo,r) - 0)1 


= 


L^(^) = 


Exgr l<^(x,0)| 


ExGr[‘^(^ + T) - (p(x, 0)]2 


ExerV5(x.O)^ 


(15a) 

(15b) 


The test case is repeated with successively refined grid 
spacing, under convective scaling (St ~ S^) and diffu¬ 
sive scaling (St S^). The exact parameterizations and 
the resulting errors are collected in Tab.|lj Fig. 11 shows 
dependency of the respective errors on the lattice res¬ 
olution. The indicated convergence rate is below first 
order, independent of the used scaling and error norm. 
This means that one cannot assert first order convergence 
with the present advection scheme in a simple translation 
test. 

With respect to the moving Laplace bubble test, the 
increased error observed in Fig. a s compared to the 
static one presented in Sec. E Al[ as well as the de¬ 


generated bubble shapes after advection (cf. Fig. |10[ ) 
suggests the following explanation. Any errors in the 
indicator function affect also the curvature computation, 
which explains the temporal oscillation of the error as the 
bubble moves relative to the grid (cf. Fig.[^. The recon¬ 
struction methods (TR and LSQR) seem more sensitive 
to the advective errors than the simpler FD scheme, and 
hence are less stable in the dynamic case. Even though 
the curvature estimates of LSQR and TR are more accu¬ 
rate and effectively reduce spurious currents in the static 
benchmarks, the combination with the low-order advec¬ 
tion scheme is problematic. 


B. Droplets on wetting boundaries 

J. Equilibrium sessile droplets 

Next, we evaluate the error at the contact line with 
solid boundaries. We change the setup of Sec. |III A] to a 
spherical cap shaped droplet, such that the initial state of 



FIG. 11: Dependency of and - errors in the 
indicator function ip on lattice resolution under 
convective and diffusive scaling (cf. Tab.|T]). Test case 
consisted of a spherical bubble of diameter d is advected 
over a distance d. The indicated order of convergence is 
below 1. 


the droplet is close to the ideal equilibrium. The equilib¬ 
rium is a spherical cap resting on the wall, with a sphere 
radius R related to the equilibrium contact angle 9eq of 
the wall hy R = h/(l — cos 0eg). Given the volume of the 
droplet V, the ideal equilibrium height of the droplet is 



(16) 


and the ideal contact line radius r (base radius of the 
spherical cap) is 


The simulated height h* and contact line radius r* have 
to be approximated from the indicator function, as 

h* = max(xz + ip('x) — 0.5), (18a) 

xg/ 


and 


1 


Xg/c 


(18b) 


where I denotes the set of interface nodes, I D Ic the set 
of contact line nodes. Hereby, f is the local approxima¬ 
tion 


f(x) = ||x-Xoll-H ((^(x) - 0.5), (18c) 

where Xq is the ideal center of the circular contact line. 
This approximation reflects that the nodes x G / are the 
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d/5x 

Sx 

St 

convective 

uo,xSt/Sx 

St 

diffusive 

uo,xSt/Sx L\ip) L^{<p) 

10 

1 

1 

0.05 

3.00% 9.38% 

1 

0.05 

3.00% 9.38% 

20 

1/2 

1/2 

0.05 

1.63% 7.19% 

1/4 

0.025 

1.38% 6.17% 

40 

1/4 

1/4 

0.05 

1.25% 7.38% 

1/16 

0.0125 

0.85% 5.39% 

80 

1/8 

1/8 

0.05 

1.17% 8.44% 

1/64 

0.00625 

0.68% 5.87% 


TABLE I: and - errors in the indicator function tp due to advection. The translation of a spherical bubble of 

diameter d over a distance d in uniform velocity field (constant velocity Ug along x-axis) is evaluated for different 
time and grid spacings St and Sx- The test case has been performed for both convective and diffusive scaling. 


centers of the corresponding interface cells. Furthermore, 
the error in the contact line is evaluated using the 
error definition, 


L^r) = 



(19) 


a. The error analysis involved the ideal contact an¬ 
gles Oeq = 30°, 45°, 60°, 90°, 120°, 135° and 150° at a 
constant resolution with Sx = 1/96 and St = 10“"^. In 
all cases, the droplet volume V was chosen equal to that 


of a hemisphere of radius R = 0.125. Fig. 12 shows the 
relative errors in height of the droplet shape obtained in 
the simulations, and the relative error in the simu¬ 
lated contact line radius. For the most extreme contact 
angles O^q, the droplet height and contact line move away 
significantly from the initial (ideal) equilibrium position, 
increasing the respective errors until the numerical equi¬ 
librium is reached. Fig. m also shows the standard de¬ 
viation in the measured contact line radius, computed 
over all contact line cells. While the errors in the con¬ 
tact line seem to be comparable in size, the higher val¬ 
ues in ST DEV (r*) obtained with the FD scheme indi¬ 
cate that the simulated contact lines are more anisotropic 
than with the LSQR scheme. 

b. A convergence study was performed for the se¬ 
lected equilibrium positions of Ogq = 60°, 90° and 120° 
by altering the grid spacing to Sx = 1/48, 1/96 and 1/144, 
using diffusive scaling for the time step. Fig. [T3| compares 
the errors obtained by the FD and the LSQR scheme 
and shows that the LSQR error is generally smaller. The 
error behavior in terms of grid dependency is somewhat 
irregular for both schemes. However, at least for the 
LSQR model the convergence rate of the error appears 
to be approximately first order. Since the construction 
of the contact points described in Sec. |II C| exploits the 
linear approximation of the reconstructed interface, the 
observed first order error is in accordance with the ex¬ 
pected behavior. 


regime (low Ohnesorge number), a spherical droplet in 
contact with an adhesive plane substrate, will start to 
spread according to the power law 




0.5 


( 20 ) 


where r is the radius of the circular contact line [53]. A 
numerical simulation of contact line dynamics requires 
the accurate modeling of a slip condition to resolve the 
stress singularities in the moving contact line. Further¬ 
more, it is in general not sufficient to work with a static 
contact angle model that imposes the equilibrium contact 
angle 6eq everywhere at the contact line |55j . In the pre¬ 
sented FSLBM, the interface representation by volume 
fractions introduces a certain amount of numerical slip 
that allows the free surface to move in the contact line 
[56] . This numerical slip is related to the grid spacing 
and does not necessarily recapture the correct physics. 
We have not introduced a dynamic contact angle model. 

We simulate the spreading of a droplet of radius i? = 10 
lattice units on a flat plate with ideal equilibrium contact 
angle O^q = 85°. The droplet is initialized as a sphere 
placed at a distance of R from the solid boundary. Due 
to a discretization effect, the simulated initial contact 
line has a positive radius r > 0, such that the adhesive 
boundary condition can be imposed in the contact line 
cells. The surface tension constant is chosen a = 4.3189 x 
10^ for a fluid of lattice reference density p = 1.0 and 
kinematic viscosity u oi vi = 3.32226 x 10“^ (first run, 
r = 0.509967) and V 2 = 1.66113 x 10“^ (second run, r = 
0.549834). Fig. 14 shows the contact line radius of the 
simulation over time. Both the FD and the LSQR model 
seem to recapture approximately Eq. (20), however, with 


a smaller exponent < 0.5. This is acceptable, considering 
the low grid resolution and the static contact line model. 
The power law obtained, r ~ seems similar to the 

one reported in m for level set-based simulations. 


IV. CONCLUSION 


2. Droplet spreading on wetting boundaries 

Because of the limitations described in Sec. mi for 
extreme contact angles, 9eq < 45° or 9f.q > 135°, we 
focus our study to a few dynamical cases. For the inertial 


Three different ways to compute the interface cur¬ 
vature from a VOF indicator function have been real¬ 
ized for comparison within a free surface LBM. A sta¬ 
tionary Laplace bubble benchmark shows that methods 
based on geometric reconstruction (TR and LSQR, in 
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FIG. 12: Sessile droplet equilibrium on solid surface for various equilibrium contact angles O^q simulated with grid 
spacing Sx = 1/96. Deviation of numerical equilibrium from ideal equilibrium for FD and LSQR scheme. 



(a) The relative deviation of the simulated droplet height h*. 



contact angle 9^q 

(b) The relative deviation of contact line radius in the 
norm, and the standard deviation of r* computed over all 
contact line cells as a measure for anisotropic artifacts. 


FIG. 13: Grid convergence study of the error in the simulated contact line radius - FD approach and LSQR 

approach in comparison. 



(a) FD scheme. 


(b) LSQR scheme. 


the present study) can reach a significantly higher accu¬ 
racy than continuum surface /orce-like approaches that 
are based on finite difference approximations (FD, in 
the present study). In accordance with previous stud¬ 
ies, reconstruction-based methods significantly reduce 
the magnitude of spurious currents. The LSQR approach 
shows a second order rate of convergence with respect to 
grid spacing, which could not be achieved with the other 
two approaches. 

For the generalized Laplace bubble test in a moving 


frame of reference, a previous study [32] reports con¬ 
vergence of errors, using a combination of higher or¬ 
der advection scheme and curvature reconstruction in a 
Navier-Stokes discretization. This behavior could not be 
reproduced with the present FSLBM. Our results indi¬ 
cate that this lack of convergence is caused by the sim¬ 
plified advection of the indicator function that does not 
take into account any geometry information. In partic¬ 
ular, the advection scheme fails to converge in a simple 
uniform advection test, indicating a lower order of accu- 


error / standard deviation 
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FIG. 14: Contact line dynamics of the spreading 
droplet; comparison of simulations based on the FD 
scheme and the LSQR scheme. The test case is repeated 
for two different viscosities, vi = 3.32226 x 10“^, and 
vi < V 2 = 1.66113 X 10“^. Both schemes tend to 
underestimate the ideal spreading law of r{t) ~ 


racy than previously reported for simple reconstruction 
based schemes (typically first or second order, with SLIC 
or PL/G-based advection in [SS]). Not surprisingly, the 
method thus fails to converge in the Laplace benchmark 
when conducted in a moving frame of reference including 
advection of the interface. This means that a conclusion 
drawn in m must be corrected: More accurate curva¬ 
ture estimation does not necessarily improve the FSLBM 
in surface tension driven flows, since (asymptotically) the 
dominant error is caused by the advection scheme. Fur¬ 
thermore, it turned out that, in the moving case, the 
curvature estimation by finite differences (FD) is less sen¬ 
sitive to errors introduced by the advection scheme than 
the reconstruction-based approaches (LSQR and TR). 

We have successfully extended the LSQR model to ad¬ 
hesive boundaries, and compared it to the FD model in 


several numerical test cases. The order of convergence 
decreases to one in the presence of adhesive boundaries. 
However, for the stationary case, the errors of the new 
approach are still significantly smaller then those of ob¬ 
tained with the non-reconstructive FD approach. In dy¬ 
namic scenarios, like contact line spreading on wetting 
surfaces, both models can recapture the basic inertial 
power law dynamics. However, similar to the bulk dy¬ 
namic case, the error situation changes. Because the 
LSQR method is more sensitive to errors stemming from 
the FSLBM advection than the FD scheme, simulations 
do not profit from the higher-order scheme. 

We conclude that a major improvement to the FSLBM 
would consist in a replacement of the advection scheme. 
VOF advection based on geometric reconstruction is sig¬ 
nificantly more complex and was, for this reason, ex¬ 
cluded from the present study. In principle, any known 
advection scheme for VOF indicator functions mi] could 
be used to replace the simplified advection scheme of the 
FSLBM [SD]. Considering that the simulation of sur¬ 
face tension according to the presented TR and LSQR 
schemes is based on a FLIC scheme to reconstruct the 
interface geometry, switching to a PLIC-based VOF ad¬ 
vection scheme would be a possible solution. A major 
benefit in accuracy can be expected. 
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